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,,_LWe study the response of a two-dimensional hexagonal packing of rigid, frictionless spherical grains due to a 
Vertically downward point force on a single grain at the top layer. We use a statistical approach, where each 
Configuration of the contact forces is equally likely. We find that the response is double-peaked, independantly 
tJfof the details of boundary conditions. The two peaks lie precisely on the downward lattice directions emanating 
dxom the point of application of the force. We examine the influence of the confining pressure applied to the 
p acking. 
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3 INTRODUCTION 



CNforce transmissions in (static) granular packings 
Khave attracted a lot of attention in recent years 



Jaeger & Nagel 1996; Bouchaud 200 3} IMueth 1998 



lair 20011 IVanel et al. 19991 |C oppersmith 1996 
edderman 19821 |Goldenberg& Goldhirsch 2002 



jeng 200 H |Geng 2003jT Granular packings are 
^assemblies of macroscopic particles that interact only 
"^yia mechanical repulsion effected through physical 
^Qntacts. Experimental and numerical studies of these 
Systems have identified two main characteristics. 
r^First, large fluctuations are found to occur in the 
Clnagnitudes of inter-grain forces, implying that the 
probability distribution of the force magnitudes is 
Bather broad (IMueth 19981 IBlair 20011) . Secondly, 
. £jhe average propagation of forces — studied via the 
Response to a single external force — is strongly 
dependent on th e underlying contact geometry 
HlVanel et al.T999l|Geng 200Tl|Geng 20031 ). 

Most of the available theoretical models capture 
either one or the other of these two aspects. The 



scalar g-model (Coppersmith 1996) reproduces the 
observed force distribution reasonably well, but yields 
diffusive propagation of forces, in conflict with exper- 
iments (Geng 2001 ). Similarly, continuum elastic and 
elasto-plastic theories ([Nedderman 1982) have been 
used by engineers for years, but they provide no in- 
formation on the distribution of force magnitudes. 
More ad-hoc "stress-only" models (Bo uchaud 2 003) 



include structural randomness, but its consequences 
on the distribution of forces are unclear. 

A simple conjecture, which provides a fundamental 
principle for the study of both fluctuations and prop- 
agation of forces, has been put forward by Edwards 
years ago (Edwa rds & Oakeshott 19891 . The idea is to 
consider all "jammed" configurations equally likely. 
A priori, there is no justification for this ergodic 
hypothesis, but its application to models of jam- 
ming and compaction has been rather successful 
(Makse & Kurchan 2002J). Its extension to forces in 
granular packings is in principle straightforward: sets 
of forces belonging to all mechanically stable con- 
figurations have equal probability. However, in an 
ensemble of stable granular packings, two levels of 
randomness are generally present (Bouchaud 2003). 
First, the force geometry clearly depends on the un- 
derlying geometrical contact network, which is differ- 
ent in different packings. Secondly, randomness in the 
values of the forces is present even in a fixed contact 
network, since the forces are not necessarily uniquely 
determined from the contact network. Instead of con- 
sidering both levels of randomness simultaneously, a 
natural first step is thus to obtain the averages for 
a fixed contact geometry, and then possibly average 
over the contact geometries. 

Such a method has been shown to produce 
single inter-grain force probability distributions in 
fixed geometry that compare well with experiments 
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(a) 

Figure 1 . Our model: N x P array of hexagonally close-packed 
rigid frictionless spherical grains in two-dimensions (drawn for 
odd N). At the top, a uniform vertical pressure p2 is applied on 
all grains, and an overload F — p2 is added on one grain. On 
the sides, a uniform pressure p\ is applied on the packing (little 
gray circles appear on interfaces where the contact forces are 
non-zero). 



( |Snoeijer 2004| ). Following a similar path, we studied 
the response of a two-dimensional hexagonal pack- 
ing of rigid, frictionless spherical grains (see Figd) 
to a pointlike external force F, and found that it is 
concentrated mainly on the two lattice directions em- 
anating from the point of application of the force 
(Ostojic & Panja 2005), in agreement with experi- 
ments ( Geng 200 Our definition of the response of 



the packing is G{i,j) = (W id ) - (W^) /F where 

W itj and W$ are the vertical force transmitted by 
the (i,j)th grain to the layer below it respectively 
with and without the external overload F. The angu- 
lar brackets denote averaging with equal probability 
over all configurations of repulsive contact forces in 
mechanical equilbrium. 

Below we summarize the method and the results 
regarding the effects of confining pressure. More pre- 
cisely, we study two cases: (i) p 2 = 0; (ii) pi=p 2 = p. 

2 CONCEPTUAL MODEL 

To start with, we describe a method for assigning the 
uniform probability measure on the ensemble £ of 
stable repulsive contact forces pertaining to a fixed 
geometrical configuration of P rigid, frictionless two- 
dimensional disks of arbitrary radii. The directions of 
the forces are fixed at each of the Q contact points, 
and one can represent any force configuration by a 
column vector F consisting of Q non-negative force 
magnitudes {F k } (with k = 1, . . . , Q) as its individual 
elements. These elements satisfy IP Newton's equa- 
tions, which can be represented as A • F = F ext . Here, 
A is a 2P x Q matrix, and F ext is a 2P-dimensional 
column vector representing the external forces. If we 
assume 2P < Q — as in the hexagonal packing we 
consider — then there is no unique solution for F. 
In that case, there exists a whole set of solutions 
that can be constructed via the three following steps: 
(1) one first identifies an orthonormal basis {F^} 



(7 = 1, . . . , d# = Q — 2P) that spans the space of 
Ker(A); (2) one then determines a unique solution 
F(o) f a ■ F(°) = F ext by requiring F<°>.FW = for 
/ = 1, . . . , die', and (3) one finally obtains all solutions 

Q-2P 

of A • F = F ext as F = F<°) + £ /, F» where 

i=i 

for I = 1 . . . Q — 2P, are real numbers. This implies 
that £ is parametrized by the //'s belonging to a set S 
obeying the non-negativity conditions for all forces, 

Q-2P 



i.e. F 



(0) 



£ /« F£ >0,VA; = 1...Q. The uniform 

i=i 

measure on £ is thus equivalent to the uniform mea- 
sure d l i = Y[ dF k 5(A ■ F— F ext ) G(F k ) = U dfi on S. 
k i 



3 ANALYSIS 
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Figure 2. Schematically shown forces on the jth grain in the 
ith layer: (a)i = 0, W^t = Vi + Sj,j F (b) i < N, (c) i = N; 
F^' j) > Vm. 



We will now apply the general method described 
above to a hexagonal packing of monodisperse, mass- 
less, rigid and frictionless disks subject to uniform 
confining pressures p\ from the sides and p 2 from the 
top (see. Fig.[l]). Our aim is to calculate the mean of 

W M = F[ id) + F 2 [i ' j) (see Fig.0 when an over- 
load F is applied on the top of the packing. 

Force balance on individual grains (see Fig. 13) can 
be expressed as 
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Note that these equations can be solved grain by 
grain for increasing i and j starting from the top 
layer. For i — 1 and j — 1, the top and side forces are 
known from the boundary conditions: W^ xt = p 2 and 
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^3 1 Pi ■ Then, if F 4 ' is chosen randomly, F 5 ' 
and Fg 1 ' 1 ^ are uniquely determined from (HJ). Simi- 

(1 2) 

larly, for the next grain i = 1 and j = 2, once F 4 ' ' is 



(X 2 1 ) i i 

fixed, F5 ' and -F 6 V1Z; are uniquely determined from 



i(1.2) 



W, 2 



p 2 , and K 



(1,2) 



F 4 . This procedure can be 



carried on identically for all the other grains for i = 1, 

except for j = P where F^°' P ^ is set by the boundary 

condition. In this way, the top forces F^'^ and F^ 
on all grains in the layer i = 2 are known, and Eq. 
® can be solved in a similar fashion to determine the 
forces for all layers successively. 

Such a sequential procedure shows that, once the 

forces F± are fixed, all the others are uniquely de- 

(i i) 

termined. In other words F 4 s for 1 < i < N, 1 < 
j < P parameterize S. The number of these pa- 
rameters is indeed dx, as it should be. Clearly, in 

(N-l,P) 

this formulation dfi = H gLF 4 (m) with F 4 {m) > 0. 

(*,i)=(i,i) 

Moreover to respect the non-negativity conditions for 
F^'^'s and F^'^'s, the parameters F±'^ must satisfy 



_ F (iJ) < F i,j 



FT 



< F (M) ; 



(2) 



implying that the set S of allowed values of F±''s is 
bounded and the uniform measure on S well-defined. 

4 g-COORDINATES AND COMPUTATIONAL 
SCHEME 



In order to evaluate (Wi ,-) 
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w^HdF, 



(k,i) 



k! 



where J\f — f s rj 4 . <iF 4 3 is the normalization con- 
stant, we define 



\/3(F 4 (ij) - if j) + F?' j) )/2] /Wi 



(3) 



where q^j is the fraction of W iy j that the (i,j)th 
grain transmits to the layer below it towards the 

left, i.e., F 5 (M) = 2g iij W i)j /V3 and F^' j) = 2(1 - 
qi,j)Wij/ \/3. Clearly, W j are the external forces ap- 
plied on the top layer. For i > 0, Wjj is a function of 
(7fc,z for k < i, since 
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The probability measure on the new coordinates gjj 
is given by the Jacobian of the variable change © 



2 NP l\dq i , j W i , j (q)/3 NP / 2 . (5) 



There is a fundamental difference between the model 
we study and the g-model (Cop persmith 1996] ): while 
in the g-model the g's corresponding to different 
grains are chosen to be uncorrelated and identically 
distributed, the uniform measure on S translates to 
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Figure 3. Simulation results forj»2 = andpi > -^F. Behavior 

of G in reduced co-ordinates x and z: (a) scaling of G(x, z) with 
system size for |x| < z/2 at two z values; (b) data collapse for 
G(x, z)\x=z/2 at three N values. 

the measure fj- Wi t j(q), which is not factorizable in 
a product of identical factors depending on a single q. 

The non-negativity conditions © for and 
Fq are automatically satisfied for < g^- < 1. 

(i i) 

However, the conditions F 4 > must be taken 
into account separately; they may further restrict the 
choice of available values of the g^j's. 

Using the g-coordinates, the (Wi t j) can be eval- 
uated numerically with a variant of the Metropo- 
lis algorithm (INewman & Barkema 1999), the usual 
Boltzmann factor being replaced by [T. . W iy j. At each 
Monte Carlo step, a grain is chosen randomly, and 
the corresponding value of q is drawn uniformally 
between and 1. The values of W^j and of the lat- 



(i i) 

eral forces F 4 are evaluated on the affected grains. 

If any F 4 becomes negative, the move is rejected, 
otherwise the new configuration is accepted with the 
usual Metropolis acceptance ratio. 

5 RESULTS 

We first consider p 2 = 0. In that case, the q^/s are al- 
lowed to take all values between and 1 if F/pi < 

The response G{i,j) then scales linearly with F 
(hence we use F = 1), and G(i,j) = outside the tri- 
angle formed by the two downward lattice directions 
emanating from (0, j ), the point of application of F. 

Our simulation results for G(i,j) within the trian- 
gle appear in Fig. |3] We find ViV that (i) the G(i,j) 
values display two single-grain-diameter-wide sym- 
metric peaks that lie precisely on the boundary of the 
triangle, (ii) the magnitudes of these peaks decrease 
with depth, and a broader central peak starts to ap- 
pear, (iii) only a central maximum for G(i,j) remains 
at the very bottom layer (i = N). 

We now define x = (j — jo)/N and (j — j + 
1/2) /iV respectively for even and odd i, and z = i/N 
(see Fig. [TJ in order to put the vertices of the tri- 
angle formed by the locations of non-zero G(i,j) 
values on (0,0), (-1/2, 1) and (1/2,1) ViV. We de- 
note by G(x, z) the response as function of these 
rescaled coordinates. The excellent data collapse in- 
dicates that the G(x,z) values for \x\ < z/2 scale 
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Figure 4. Simulation results for G(i,j), for pi = P2 = 1: (a) 
G(i,j) for F = 10 as function of j — jo> at f° ur different values 
of i for various system sizes; (b) G(i,j) as function of i for jo — 
j = i, for F/p = 5, 10 and 15 and three different system sizes. 

with the inverse system size [Fig. Ha); we how- 
ever show only two z values], while the G(x, z) val- 
ues for \x\ = z/2 lie on the same curve for all sys- 
tem sizes [Fig. Htb)]. The data suggest that in the 
thermodynamic limit iV — > oo, the response G(x,z) 
scales ~ 1/N for \x\ < z/2, but reaches a non-zero 
limiting value on \x\ = z/2\fz < 1. We thus ex- 

p6Ct i^o G ^' Z )l|^l=^/ 2 > G ( x '*)Im<*/2 V2: < ^o 1 " 
equivalently, a double-peaked response at all depths 
z < 1 in the thermodynamic limit. 

For F/pi > & and p 2 = 0, as the value of p\ de- 
creases, the qi/s get restricted to narrower ranges 
within (0, 1) with increasing values of F/pi, and con- 
sequently, an increasing amount of vertically down- 
ward force carried by the grains is transferred to the 
boundary of the triangle. In the limit F/pi — > oo, the 
non-negativity conditions of all the forces make all 
F^'^ and F^'^ vanish, and the packing effectively 
becomes rectangular. In fact, the same behaviour also 
occurs for any value of p 2 ^ 0. 

For p 2 7^ however, experimentally the most rele- 
vant case is p% = p 2 = p, since in practice a confining 
pressure has to be applied on the packing. In this case, 
the response clearly depends on the dimensionless pa- 
rameter F/p. The simulation results for iV = 20,30 
and 40 are presented in Fig.EJ 

Fig. |4](a) shows the response for F/p = 10 as func- 
tion of j — jo for various values of i. For small i, 
the response displays two single-grain-diameter-wide 
symmetric peaks along the lattice directions emanat- 
ing from the point of application of the overload. As i 
increases, the magnitude of the peaks decreases, and 
the response for large i becomes essentially flat. The 
dependance of the response on the system size in this 
case however is very different from the case p 2 = 0. 
Rather then being self-similar, the response is inde- 
pendant of N: G(i,j) at a given depth i is the same 
for all systems with N > i layers. 

The values of G(i,j) on the lattice direction i = 
jo — j as function of i are plotted in Fig.Etb), for N = 
20, 30 and 40, and F/p = 5, 10 and 15. For increasing 
F/p, the values G(i,j — i) increase for fixed i and 
their decay with increasing i is slower, implying that 
the peaks are more pronounced and propagate deeper 



in the packing. 

6 CONCLUSIONS 

In summary, we find that assigning equal proba- 
bility to all mechanically stable force configura- 
tions for rigid, frictionless spherical grains in a 
two-dimensional hexagonally close-packed geometry 
yields a double-peaked response independantly of the 
details of the boundary condition. The peaks are sin- 
gle grain diameter wide, and they lie on the two down- 
ward lattice directions emanating from the point of 
application of F. In the case of zero top pressure, the 
response exhibits self-similarity, but when the pack- 
ing is confined by uniform pressure, the peaks pene- 
trate the packing deeper with larger F. Such a simple 
model is in good qualitative agreement with experi- 
ments in a reasonably robust manner. Whether grains 
with friction (Breton 2002) produce broadening of the 
peaks or not remains to be investigated. 
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